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Abstract 

In order to model pressure and viscous terms in the equation for the Lagrangian dynamics of the velocity gradient tensor in 
turbulent flows, Chevillard & Meneveau {Phys. Rev. Lett. 97, 174501, 2006) introduced the Recent Fluid Deformation closure. Using 
matrix exponentials, the closure allows to overcome the unphysical finite-time blow-up of the well-known Restricted Euler model. 
However, it also requires the specification of a decorrelation time scale of the velocity gradient along the Lagrangian evolution, and 
when the latter is chosen too short (or, equivalently, the Reynolds number is too high), the model leads to unphysical statistics. In 
the present paper, we explore the limitations of this closure by means of numerical experiments and analytical considerations. We 
also study the possible effects of using time-correlated stochastic forcing instead of the previously employed white-noise forcing. 
Numerical experiments show that reducing the correlation time scale specified in the closure and in the forcing does not lead to a 
commensurate reduction of the autocorrelation time scale of the predicted evolution of the velocity gradient tensor. This observed 
inconsistency could explain the unrealistic predictions at increasing Reynolds numbers. We perform a series expansion of the matrix 
exponentials in powers of the decorrelation time scale, and we compare the full original model with a linearized version. The latter 
is not able to extend the limits of applicability of the former but allows the model to be cast in terms of a damping term whose sign 
gives additional information about the stability of the model as function of the second invariant of the velocity gradient tensor. 

Key words: turbulent flows; turbulence simulation and modeling; isotropic turbulence; homogeneous turbulence 
PACS: 47.27.-i, 47.27.E-, 47.27.Gs 



1. Introduction 

The study of the velocity gradient tensor Aij = dui / dxj , 
where u{x, t) is the velocity field, has great significance in 
determining the small-scale structure of turbulent flows. 
Among others, it allows to identify the areas of the flow 
in which either strain-rate (i.e. deformation) or vorticity 
(i.e. rotation) prevails. It contains geometric information 
about the orientation of vorticity and strain-rate eigen- 
vectors, and its higher-order moments quantify the level 
of intermittency of small-scale turbulence. Here we focus 
on three-dimensional, incompressible turbulent flows, i.e. 
d ■ u = An = 0. 

The Lagrangian interpretation of the evolution equation 
for A, derived from the (Eulcrian) Navier-Stokes equation, 
leads to a closure problem due to the fact that the pres- 
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sure Hessian and viscous term are not directly expressed in 
terms of the local value of A. The Restricted Euler (RE) 
model |l|2j is the classical and crudest closure scheme, as 
it completely neglects the viscous term and the deviatoric 
part of the pressure Hessian, assuming the latter to be 
isotropic. The RE closure captures several qualitative fea- 
tures of real turbulence, such as the alignment of the vor- 
ticity vector with the direction corresponding to the in- 
termediate eigenvalue of the strain matrix, and the statis- 
tical prevalence of axisymmetric expansion (situations in 
which a fluid element, initially spherical, become pancake- 
like rather than cigar-like). It allows the analytical investi- 
gation of the Lagrangian evolution of material volume el- 
ements and of their deformation tensor |3j. Moreover, it 
also represents the starting point to deduce the so-called 
Advected Delta- Vee (ADV) system for velocity increments 
|4|5j. However, it presents a major drawback, namely the 
appearance of an unphysical finite-time singularity. 

With the aim of overcoming this problem, the so-called 
Recent Fluid Deformation (RED) closure was introduced 
|6|7|8j . It provides a model for the pressure Hessian and 



viscous term, which are suitably parameterized, i.e. ex- 
pressed as functions of the local velocity gradient tensor. 
Using matrix exponentials, it takes into account both the 
geometry and the dynamics of the recent history of the 
deformation of a fluid element, and requires the specifi- 
cation of a decorrelation time scale r, of the order of the 
Kolmogorov time scale. It was shown that this refined 
parameterization leads to stationary statistics, removing 
the above-mentioned finite-time singularity. Moreover, it 
also captures some other relevant properties observed in 
real turbulence, such as: the non-Gaussianity of longitu- 
dinal and transverse velocity gradients, the characteristic 
teardrop shape of the probability density function (PDF) 
in the R-Q plane, the quasi-lognormality of pseudodissipa- 
tion, and, more importantly for our purposes, the correct 
scaling with the Reynolds number (Re ~ t"^) at small 
to intermediate values of Re. However, when t is chosen 
too short (or, equivalently. Re is too high), unphysical 
statistics are observed in the model [617] . 

In the present paper, we investigate the limitations of the 
RFD closure at increasing Reynolds numbers (or decreas- 
ing imposed correlation time scale, r) by means of numer- 
ical computations. In particular, the predicted Lagrangian 
time-correlation structure of the model is studied. We also 
propose a simplified, linearized version of the model, which 
turns out to have the same (or even narrower) ranges of 
applicability, but allows us to proceed analytically (in the 
spirit of [S]) and to draw some additional conclusions. 

The paper is organized as follows. In scction[2]we review 
the RE model and RFD closure. The equation deduced from 
RFD is then integrated numerically in time, and the results 
are shown in scction[3l In section[4]we perform an analytical 
manipulation of the matrix exponential, expanding it in a 
power series in r, and we show the resulting equations for 
A at the different orders in t. In section[5]we focus on the t- 
linear approximation and we deduce a restricted dynamical 
system. Conclusions and perspectives follow in section [6l 
Appendix [Al shows some details about the random forcing 
term we use to obtain stationary statistics for the RFD 
closure. 

2. Review of RE model and RFD closure 

Starting from the Navier-Stokes equation, 
dtUi + UkdkUi = -dip + i^d'^Ui 

{p is the pressure divided by density and u the kinematic 
viscosity), it is sufficient to take the gradient {dj) to obtain 
the evolution equation for the velocity gradient tensor: 

dtA,j + UkdkAij = -AikAkj - dtdjp + lyd'^A^j . (1) 

The incompressibility constraint makes A a traceless ten- 
sor, tr j4 — 0, thus reducing its number of independent 
components from 9 to 8. 

Equation ([1]), interpreted in Eulerian sense, tells us little 
more than the original Navier-Stokes formulation. The key 



point consists in reinterpreting it in a Lagrangian sense, by 
considering its left-hand side as a material derivative fol- 
lowing a fluid particle. However, one immediately faces a 
particular closure problem, because the last two term at 
its right-hand side (RHS) are not known in terms of A at 
the same point and time. The isotropic part of the pressure 
Hessian can easily be rewritten by using the Poisson equa- 
tion (obtained from ^ by taking its trace), which allows 
to express the pressure Laplacian as a quadratic term in A: 
d^p = —trA^. Thus, one obtains: 



where 



AikAkj ~ AkiAik-^ I +H. 



(2) 



Hi 



d,djp - 



d^p^-^ 
3 



vd'^A, 



is the (traceless) unclosed term for which a suitable model 
is required in order to express it in terms of A and other 
known quantities. Notice that, if one wants to investigate 
not the actual velocity field itself, but rather a coarse- 
grained version of it (e.g. in order to study velocity incre- 
ments, as we will show at the end of section [SJ in the spirit 
of a large-eddy simulation (LES)), equation remains 
valid, provided one includes the appropriate subgrid terms 
into [TUTTT] . 

The RE model takes H — 0, i.e. it completely neglects 
the effect of viscosity and the anisotropic action of the (Eu- 
lerian) pressure Hessian, which are known to play an im- 
portant role in the regularization and evolution of the tur- 
bulent dynamics. It is therefore not surprising that the RE 
model 



dtA 



trA^ 



leads to a finite time singularity. 

The RFD closure overcomes the blow-up by introducing 
four approximations and two time scales, namely a typical 
decorrelation (Kolmogorov) time r and an integral (large- 
eddy turnover) time T. The closure is based on the mapping 
between the Eulerian position x (at time t) and an initial 
Lagrangian coordinate (or label) X (at some earlier time 
t — t). This map is invertible because of incompressibility, 
as the Jacobian has unit modulus. We briefly recall the 
four approximations used to obtain the closure. For more 
details, see [STTlg] . 

As a first step, when writing second derivatives by means 
of the chain differentiation rule, one neglects the spatial 
variations of the Jacobian: 



dp _ dXm dp 
dxi dxi dXr,t 

dAjj _ dX,n dAjj 

dxk 



d^p dX„, dXn d^p 



dxidxj dxi dxj dXmdXn 
i/d^Aij dXm dXn vd'^Aij 
dxk dXra dxkdxk dxk dxk dXrndX„ 



Secondly, the Lagrangian pressure Hessian is assumed 
isotropic (this is physically more meaningful than the cor- 



responding approximation in RE about the isotropy of the 
Eulerian Hessian): 



oc S„ 



dXjndXn 

(the proportionality prefactor can be found through the 
Poisson equation). Moreover, one models the viscous term 
via an isotropic linear damping with the integral time scale 
T: 



dXmdXn 



4 ■ A 

^ 3~ 



Finally, the crucial step is represented by the parameteri- 
zation of the well-known Cauchy-Green tensor, 

with an "on-off approximation" : 

C ~ e^^e^^' , (3) 

i.e., the actual slow decorrelation of C along the Lagrangian 
trajectory is replaced by a perfect correlation only during 
the Kolmogorov time scale r and by a complete decorre- 
lation for longer times: in this way, the cumbersome, full 
time-ordered exponential reduces to a much simpler ma- 
trix exponential (here, the superscript ^ denotes the trans- 
posed matrix) . The Reynolds number (based on the integral 
scale) corresponds to Re — (t/T)~^, or in other words the 
Taylor-scale Reynolds number behaves as Re^ ~ r~^. From 
([2]), using these approximations, and including a stochastic 
forcing term F to enable statistically stationary dynamics, 
the resulting system is: 



dtA = -A^ + C 



,_i trA^ 



trC-i 



- A 



trC 



-1 



3r 



(4) 



with C given by Notice that the RE model corresponds 
to r = 0, such that C = I, with = 0. In what follows, we 
will thus try to understand what happens when t becomes 
smaller and smaller, a critical aspect of the RFD closure. 

Among related approaches to remove the blow-up behav- 
ior of RE, we can mention the tetrad model |12I13I14I15| . 
the Cauchy-Green models I16ll7j . and a multi-scale cou- 
pling approach in the spirit of shell models [18] . The RFD 
closure is attractive because of its relative simplicity: it is a 
system of only 8 independent stochastic ordinary differen- 
tial equations, based on a physically inspired connection to 
the pressure Hessian. However, unlike the shell model which 
allows to be extended to high Reynolds numbers by addi- 
tion of shells, as already mentioned the RFD closure pro- 
duces unrealistic predictions at increasing Reynolds num- 
bers. In } 6|7j . T ~ 0.06 was the smallest value for which re- 
alistic results could be obtained. The analysis that follows 
explores the behavior at smaller r in detail, in order to pro- 
vide a basis for possible future improvements of the model. 



3. Numerical investigation of the dynamics 

In this section we present results based on numerical time 
integration of (|4]) , including the full matrix-exponential ex- 
pression for C. In terms of documenting results from the 
simulations, one can follow the evolution of the 9 (8 of which 
independent) components of A individually, e.g. in order to 
investigate their time correlations. One can also represent 
parts of the relevant dynamics in terms of the two invariants 
R = -tr A3/3 and Q = -tTA'^/2. In the R-Q plane the 
zero-discriminant curve 27R^ + 4:Q^ = is usually called 
the Vieillefosse line (denoted in the plots as Vline): the RE 
model is known to diverge on its right-bottom part. 

As numerical procedure we adopt a standard fourth- 
order Runge-Kutta scheme for the time integration, with 
time step At small enough (always At < 10~'^t, i.e. at least 
one thousandth of the shortest relevant time scale). As in 
|6|7j , we perform a non-dimensionalization of the equations 
with T, thus the only remaining significant parameter is r 
(- Re-^/2). We run the code for a total time > 10'^ (i.e. 
at least one thousand large-eddy turnover times), in order 
to accumulate well-converged statistics. As a first test, we 
integrated the equations starting from several initial con- 
ditions without forcing {F = 0), in order to check that the 
dynamics properly decay, converging to the origin. 

In order to obtain stationary statistics, the tensorial forc- 
ing term F is introduced on the RHS of ^ . Following argu- 
ments given in |6|7| , this term mimics the effects of neigh- 
boring small-scale eddies that affect the velocity gradient in 
its Lagrangian evolution. We performed some preliminary 
tests by using a white-in-time random signal (following [6] ) , 
but then we chose to focus on a Gaussian, smooth noise, 
with a finite time correlation 9. This time scale will be cho- 
sen of the order of r itself, because this is the time scale at 
which the interactions with the other eddies are believed to 
take place. A comparison between the two kinds of noise, 
and thus on the influence of the time correlation, will be 
made in subsection 13. 2) where various values of 9 (namely 
9 = T, t/2, t/5, t/10) will be considered. In appendix [Al 
we discuss some details of the forcing. 

From ([3]) it appears that the tensor C will become im- 
portant once A should have order of magnitude ~ r^^. 
The quadratic terms on the RHS of ^ then suggest that A 
should vary on a characteristic time scale of order r. Thus, 
as a first guess for the scaling of the amplitude of the noise, 
it was chosen of the order of r"-^ so that it is comparable to 
the magnitude of A^. However, for r ~ 10~^ and smaller, 
this turned out to be a too strong a forcing amplitude: the 
forcing overwhelmed the quadratic self-interaction term, 
resulting in nearly Gaussian statistics for A. Therefore, a 
more careful procedure must be followed to prescribe the 
amplitude of the forcing. We begin by taking the value r = 
10~^ as a baseline case (this value was the baseline value 
also for the prior studies [6|7|8j ). For this reference case, in 
[6 8] an amplitude of the forcing F = 0{1) was used, along 
with a compensation with the integration time step, due 
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Fig. 1. Response behavior of multiplied by r, as a function 

of the typical noise amplitude J^, for three different values of r. The 
two (approximate) intercepts with the horizontal line passing through 
the point representative of the reference case r = 10~^ (for which 
one only value is clearly plotted, corresponding to = 1) represent 
the correct values of which will be used in what follows for the 
cases T = 10"'^ and = 10~^, i.e. J^' and JF" respectively. The arrows 
have been introduced to help the understanding of the procedure. 
For every point (full squares and circles), the corresponding couple 
of empty symbols (above and below it, in some cases overshadowed) 
indicates the minimum and maximum values obtained by subdividing 
the whole numerical run into 10 equal sub-intervals, and is thus 
an indication of the statistical uncertainty appearing in the ranges 
where the increase in amplitudes is very steep. 



Fig. 2. Joint PDF of R-Q obtained from the RFD closure, using 
r = 10~^ and forcing with time correlation 9 = t and typical 
amplitude = I. Note that P{R, Q) is plotted as a function of 
scaled values t^R and t-^Q, instead of R and Q. 
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to the white-noise character of the forcing used there. Let 
us denote the actual amphtude of the forcing Equation 
(jA.l[) in the appendix shows how the forcing is generated 
as function of J-. Since we wish to compare our results with 
the ones of |6l8j , we also adopt T = 1 for this reference case. 
As shown in appendix [Aj J- is used in (jA.ip as amplitude 
of the diffusion term. Note that Fl'^"" = T. 

As a next step, wc consider increasing Re, and we know 
that the strength of the forcing (i.e. of the action of neigh- 
boring eddies) should increase, but at a rate smaller than 
^ T~^, as already explained. An important criterion that 
the dynamics should follow when increasing Re is that, for 
energy dissipation to remain constant at increasing Re, one 

should observe that = ^J{A^ ~ r"-^, when r is 

interpreted as the Kolmogorov time scale [B]. Hence, the 
forcing may be adjusted to generate approximately such 
scaling. In figure [T] we plot, in logarithmic scale, the re- 
sponse behavior of multiplied by r, as a function 
of T , as obtained from a large number of simulations. We 
observe that t||A||'''°^ increases monotonically until a cer- 
tain forcing strength, and then a sudden increase in ampli- 
tude is observed. Beyond this transition, the increase con- 
tinues less dramatically. Following the energy-dissipation- 
based criterion (of aiming to reproduce t||A||''™*' ~ 0.08), 
the values J-' = 4.5 and T" ~ 14.35, obtained by finding 
the approximate intercepts with the horizontal line drawn 
through the reference-case point, represent the correct am- 
plitudes to be used for t = 10~^ and = 10"'^, respectively. 
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Fig. 3. Same as in figure |2] but for t = 10^"^ and T = 4.5. 

3.1. Baseline results 

The baseline case of r = 10^^ is computed, using a Gaus- 
sian random forcing with time correlation 6 = t and with 
forcing parameter — I. The model yields the probabil- 
ity density function for Q and R shown in figure [2l which 
agrees well with [6 . It is worth mentioning that, differently 
from here we multiply Q and R by (the appropriate 
power of) r, rather than normalizing them with the trace 
of [iA + A'')/2]\ 

Next, we consider the case of r — 10^^ (with T — T' ^ such 
that becomes 10 times bigger than in the previous 

case) and show the results in figure [3l The shape of the 
PDF is reasonable, but one should notice that the values of 
R and Q do not scale with factors lO'^ and 10^ with respect 
to figure [51 as one would expect from a simple dimensional 
analysis. This is a clear indication of the presence of a very 
intermittent behavior, as we will show in subsection [321 
This observation is confirmed further when plotting the 
case r — 10~^ (figure[5]), for which the picture is very sim- 
ilar to the latter case, in the sense that the rescaling of R 
and Q is much smaller than expected. 

One could therefore wonder what happens when the am- 
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Fie. 4. Same as in figure [1 but for t = 10"^ and = 14.35. 
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Fig. 5. Same as in figure[3l but for = 31.6. 

plitude of the forcing is chosen in a way such as to give re- 
sults whose orders of magnitude of Q and R are consistent 
with the expected Q ^ r^^ and R ~ scahng. In other 
words, for r = 10~^ we momentarily impose a different 
value J- = 31.6 T') such that, passing from figure[2]to 
fi-gureEl a rescaling of r by one order of magnitude (10^^) 
brings about an approximate rescaling of the quadratic and 
cubic quantities, Q and R, by factors 10^ and 10"^ respec- 
tively. The results are still reasonable in terms of the shape 
of the PDF. 

However, at t — 10^^ (figure[6]), when using the same crite- 
rion (which implies T = 1000, ^ J-") the model yields the 
results shown in figure [Sj The resulting PDF displays an 
unrealistic diamond-like shape skewed towards the second 
and fourth quadrants. Further modifying the amplitude of 
the noise does not bring about better results for this last 
case. 

Besides testing other forcing amplitudes, we have also 
tried other possible forms for the forcing, including various 
multiplicative noises. Among them, we can mention: dtA = 
... + (F-Af,dtA = ... + F||A||2,dtA = ...+F(1 + ||A||2) 
and dtA = . . . + F{1 + Q), where the superscript " identifies 
the deviatoric part of the tensor (i.e. the tensor minus its 
trace times //3). However, none of them lead to physically 
meaningful results for small r. 

Analogously, following [7j , we tried to impose not a fixed. 



0.01 



0.005 



-0.005 



-0.01 



( 


1e-10 
1e-11 
1e-12 

., 1e-13 

Viine - 







-0.001 



-0.0005 




r'R 



0.0005 



0.001 



Fig. 6. Same as in figure U but for JF = 1000. 

overall dissipative time scale, but rather a variable, local 
one, T{t) oc [tr (A + A-'^)^]"^/^. However, also this approach 
did not allow us to obtain better, or more physical, results 
for cases corresponding to higher Reynolds numbers. 



3.2. Time correlation structure of velocity gradients 

In this section, we document time series, autocorrelation 
functions and time scales of velocity gradient tensor ele- 
ments as predicted by the model. Several options of the 
forcing term are considered. Specifically, we vary the time 
correlation 6 of the forcing. Five cases are considered: the 
baseline case with 9 = t itself, the white-noise case (i.e. 
6 = 0), and three more finite-correlated instances with in- 
termediate time correlations of 6* = r/2, t/5 and r/lO. 
Each case needs a modified amplitude J-^ in order to try 
to maintain the amplitude of the velocity gradient tensor 
comparable, as 9 is varied. Inspired by the behavior of a 
forced linear system, we use the rescaling in which is re- 
placed by J-\0=t\/t/9 for the finite-correlated forcings and 
by T\0=r\/T/ At for the white-noise case. It is observed 
that this rescaling indeed maintains the typical amplitudes 
of A unchanged even as 9 is varied. 

In figure [71 we plot the signal, as a function of time, of 
one longitudinal (^ii) and one transverse {A12) component 
of the velocity gradient tensor, at r = 10^^. The same 
analysis is repeated in figures [HHH at = 10~^ and 10~^ 
respectively, only for the longitudinal component. In all 
cases, there seem to be essentially no differences between 
the five evolutions plotted, which correspond to different 
time correlations 9 of the forcing. Note from the inserts 
(corresponding to longer, full numerical runs) that the cases 
for T = 10~^ and 10~^ display some very infrequent but 
violent fluctuations, which are absent for the baseline case 
T = 10^^. These fluctuations are responsible for the sudden 
increase in amplitudes seen in figure [TJ A detailed study of 
the dynamics during such intermittent excursions is left for 
a future study. 

Next, the time correlations of A are quantified, for sim- 
ulations using the same five forcing cases. In figure [TOl for 
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Fig. 7. Time signal of one longitudinal (top panel) and one transverse 
(bottom panel) component of the simulated velocity gradient tensor, 
multiplied by r, at t = 10~^ and T = 1. Five lines are shown, 
corresponding to different time correlations of the forcing. Note that 
the plotted time window has been chosen as 100 (i.e. lOOT in units), 
so that r corresponds to one thousandth of the total horizontal 
extension. 
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Fig. 8. Same as in the top panel of figure [7] but with r = 10~^ and 
T = T' . Note that r corresponds to one tenth of thousandth of the 
total horizontal extension. The insert shows the results of the full, 
longer run, whose first part is reproduced in the outer plot. 



both one longitudinal and one transverse component, we 
plot the autocorrelation function 
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Fig. 9. Same as in figurc|8]but with r = 10"^ and = JF". Note that 
T corresponds to one hundredth of thousandth of the total horizontal 
extension. 

(no summation implied), where the temporal average has 
been performed on our entire numerical run. Figure 1111 
shows the results for r = 10~^ and 10"'^ as well for long 
iterations (only the three cases 9 = t are plotted). A strik- 
ing observation is that the measured correlation time scale 
from the autocorrelation function is quite large, and in both 
cases reaches 0.2 only at a time scale s ^ T {= 1), signif- 
icantly larger than r. More surprisingly, when decreasing 
the imposed time scale r (which corresponds to an assumed 
decorrelation time in the RFD closure), the results from 
the simulations certainly do not show a concomitant de- 
crease in the predicted autocorrelation function. Also for 
T = 10~^ and 10~^, we observe that the autocorrelation 
only approaches zero for s ~ T and not r. 

4. Power-series expansion of the matrix 
exponential 

One major effect of the RFD closure for the pressure Hes- 
sian is to oppose the generation of finite-time singularities 
that occur in the RE equation. Hence, at least in parts of the 
R-Q plane, the closure acts as a regularization or damping 
effect. Detailed comparisons with direct numerical simula- 
tions (DNS) in [8] confirm this effect using conditional av- 
eraging analysis. In order to study this issue by inspecting 
the closure terms directly, it is convenient to introduce a 
simplification of the matrix exponential using expansions. 
The definition of the matrix exponential is used to obtain 
a power-series expansion: 

e-^ = y i-(TA)" = I + tA + -t^A^ + -r^A^ + ... . 
^-^ n\ 2 6 

Replacing in ([3]) and (g]), one obtains: 
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Fig. 10. Time autocorrelation of one longitudinal, aii(s) (top panel), 
and one transverse, ai2{s) (bottom panel), component of the velocity 
gradient tensor, at t = 10~^ and for five different values of 9. Notice 
that, in the linear scale (outer plot), t corresponds to one hundredth 
of the total horizontal extension lOT. In the insert, the axis of the 
abscissae is represented in logarithmic scale, and the vertical line 
indicates s = t. 
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Fig. 11. Time autocorrelation of one longitudinal component of the 
velocity gradient tensor, aii(s), for three different values of t. In the 
insert, the axis of the abscissae is represented in logarithmic scale. 
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It is also interesting to recast ([5]) by separating it into 
its symmetric and antisymmetric components, S = {A + 
A^)/2 and n = A - 5. We get: 

dtS = - + n^Y - -Q [2{s'^ f + sn- ns] 
[iisy - 2nsn + sfi^ + fi^s 

s r 4 4,4, 



dtn 



4 c. 



(6) 



(7) 



Here, Qs = -trS'^/2 and Rs = are the second 

and third invariants of S (its first invariant, i.e. its trace, is 
identically zero). Notice, from ([6]) and ([7|), that symmetric 
initial conditions (f2|t=o = 0) only evolve into symmetric 
dynamics (pure strain), while this is not true for the an- 
tisymmetric counterpart (pure rotation), because starting 
from S\t=Q = a strain rate can develop. It is also worth 
mentioning that the symmetric pressure Hessian has no di- 
rect influence on vorticity: this is why, in ([7|), the time scale 
r only appears in the dissipative terms (which lack the or- 
der r^), while no contribution from the pressure Hessian is 
present. On the other hand, the latter is the origin of the 
0(t) term in the last line of which is not derived from 
viscosity (indeed it does not depend on T) but has been 
written in this way for the sake of simplicity. 

The simplification induced by the projection on S and CI 
allows now to identify a linear term, multiplied by the term 
inside the curly parentheses. The sign of this term can then 
be investigated: we will later return to this issue in section 
[21 in the framework of a simplified dynamics, obtained by 
discarding the 0{t^) (or higher) terms in 

4.1. Unforced dynamics 

Next, we wish to compare the numerical results of the 
full matrix-exponential approach (to be labeled as Mexp) 
with the ones obtained by integrating ([5]), in which one can 
retain e.g. up to the first, second or third order in r (la- 
belled, respectively, as lin, quad, cub). 
In figures [T2HT3] we show the unforced evolution for the full 
matrix-exponential model and for its linear, quadratic and 
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Fig. 12. Evolution of _R and Q (multiplied by the appropriate power of 
r) for the unforced dynamics using the RFD closure with r = 10~^, 
imposing the full matrix-exponential model or its linear, quadratic 
and cubic approximations. The four arbitrary initial conditions have 
been chosen as: (a) Aij\t=<^ = 3i — 2j, such that Q\t=Q = 18 and 
K|t=0 = 36; (b) Aij|t=o = — 3i + 2jr, such that Q|t=o = 18 and 



R\t=Q = -36; 
and A32|t=o = 



(c) Aij\t=Q = 3i - 2j (except for Aailt^o = -5 
-7) such that Q|t=o = -18 and R\t=o = -36; (d) 
Aij\t=o = 3i ~ 2j (except for yl3i|t=o = -5 and A32|t=o = -1) 
such that Q\t=o = — 18 and R\t=o = 36 (in all cases we have then 
redefined Axi\t=o = --4ii|t=o - A22\t=o)- 
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Fig. 13. Same as in figure [T2l using t = 10 ^, but not including case 
(d), which has not been reproduced to avoid crowding in the figure. 

cubic approximations, starting from a few representative 
initial conditions in each of the quadrants of the R-Q plane. 
Figure [T2l shows the resulting time evolution for r = 10~^, 
for which all these unforced cases converge and decay to- 
wards the origin. For r = 10~^ (figure [T5|) similar results 
can be obtained, although the weaker restitution term pro- 
vided by the lower value of r allows some excursions along 
the bottom-right Viellefosse tail before the decay towards 
the origin. These trends also hold at r = 10^'^ (not shown), 
also starting from different quadrants, but then numerical 
instabilities occur if r is further reduced, and these can only 
be avoided if the integration time step is reduced drasti- 
cally. 

The important observation is that, in all of the previous 
cases, the linear, quadratic and cubic expansions are similar 
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Fig. 14. Same as in figure [2] with both the full matrix-exponential 
model and its linear approximation shown. Proceeding from the 
interior outwards, the isolines are for densities of 10, 1, 0.1 and 0.01 
respectively. 

to one another, and provide a good approximation of the 
full matrix exponential. Therefore, in what follows, we will 
only focus on the comparison between the linearized model 
and the full, original one, not focusing on the O(t^) and 
O(t^) approximations. 

4.2. Forced dynamics 

Let us now introduce again our additive forcing term, 
with a finite correlation time 9 = t (according to the frame- 
work described in section[3]and in appendix lX)) . also on the 
RHS of ([5]) only including the first order in r. 
Kl T — 10^^, and T = \, the linear approach is basically 
equivalent to the (already analyzed) matrix-exponential 
one in terms of resulting PDF's, as shown convincingly in 
figure [HI 

Reducing r, the simulations with the linear approxima- 
tion based on the energy-dissipation criterion {T — 4.5 for 
T = 10^^ and T — 14.35 for r — 10""^) are again basi- 
cally equivalent to the ones from the full matrix exponen- 
tial shown in figures[3Hll However, if one tries to reproduce 
the dimensional-scaling-based results of figures [5H6] {T = 
31.6 for T = 10"2 and T = 1000 for r = lO^^)^ simulations 
with the linear approximation begin to suffer from numer- 
ical instabilities, regardless of how small the time step is 
chosen. It is therefore not possible to perform detailed nu- 
merical comparisons for significantly lower values of r. In 
order to proceed one step further using analytical tools, in 
the next section, the evolution equation for the complete 
set of relevant invariants is studied. 

5. Analytical study of the invariants in the 
"linearized" equations 

Focusing on the "linearized" model, derived from ((5]) by 
retaining only up to the first-order terms in r, allows us 
to proceed one step further analytically. Following [9J, we 
derive a dynamical system for the five tensor invariants 



Q, R, Qs, Rs (already introduced in section 2]), and V"^, 
the latter defined from the strain-vorticity alignment V = 



I S • u; I , where uji 



-Cijk^jk is the vorticity vector. Notice 



that Qs < < hy definition, moreover it can easily be 
shown that Qs < — (27 i?s^/4)^/'^, so that the dynamics is 
confined below the zero-discriminant curve in the Rs-Qs 
plane [S] . It is worth mentioning that it is not necessary to 
investigate the antisymmetric counterparts of Qs and Rs, 
because by definition Qfi = Q — Qs and Rfi — 9J . 

Starting from the equation for Aij (without forcing) and 
using appropriate contractions and the Cayley-Hamilton 
theorem, it is possible to obtain the following system of five 
equations: 



-3R-2T- 



dtQ 

dtQs = -R-2Rs-2T 

dtR 
dtRs 
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Notice that, in order to keep track of the various terms, we 
have not non-dimensionalized the equations and kept the 
term T . System ([8|) shows linear and quadratic couplings 
among the different variables. What is most interesting is 
the fact that the linear "damping" terms, which have been 
factored at the end of each line for the sake of clarity, in- 
clude terms at 0{t). It is of particular interest to study 
the sign of the square braces: when they (or at least one 
of them) becomes negative, divergences and blow-ups are 
likely to occur. However, in all equations except the one for 
Qs, these terms are positive definite (in fact larger than 
1) since Qs is never positive by definition. The possibility 
of negative damping is possible in the equation for Qs it- 
self, in which the damping term may become negative when 
Q > 3/{ATt), that is in highly rotating regions when Q is 
large and positive. In fact, the rate of phase-space contrac- 
tion (—V • z, where z is the vector formed by the RHS of 
©) is 14T-1 - 40Qst/3 - SQut/S, i.e. for overaU grow- 
ing solutions one would need Qn > 5|(5s| + 21/(4Tt), very 
large rates of rotation. Clearly, however, along the excur- 
sions of the system along the Vieillefosse tail (the right- 
bottom part of the R-Q, where the RE model is known to 
diverge) when Q is negative, the term provides damping 
and thus protects against divergencies along this direction. 

To obtain a system that generates stationary PDF's it 
is difficult to introduce random noise properly on the RHS 
of ([8]). The difficulty stems from the inherent bounds and 
internal consistency conditions among the tensor invariants 
(e.g. the limits on Qs and V'^ mentioned above). Simply 
adding random noise to the equation system[5]quickly leads 
to Qs > or < 0, etc.. Introducing a simple, additive 
random forcing on the RHS of the evolution equation for 
A, dtA — . . . + F, leads to unclosed terms on the RHS 



of ([H]). Therefore, the analytical study of the invariants is 
limited to the results presented above. 

Another set of interesting variables can be constructed 
by considering the tensor invariants of A combined 
with a material line element r being convected by the 
flow. Specifically, following [4l5j . one can introduce 
Su = iAijriVj/r'^ (which can have either sign) and Sv = 
\£{Sij — TiTj /r'^)Ajkrk/r\ (the magnitude of the transverse 
velocity vector Svi). These two variables represent the lon- 
gitudinal and transverse velocity increments at a fixed scale 
i, i.e. between two fixed points with constant separation 
(rather than between two material points or fluid particles 
whose relative distance r{t) would change with time). The 
derivation of a Lagrangian time evolution equation for Su 
and Sv based on the equations for A and r leads to the 
ADV system introduced and studied in 4:5J. In this prior 
work, it was shown that the neglect of viscous and much 
of the pressure Hessian term leads to unphysical behavior 
of the dynamics. Hence, it is of interest to consider the 
implications of the RFD closure and its expansions in the 
context of the ADV system. 

Repeating the relevant derivations [4 5J but including the 
RFD closure expanded to first order in the equation for A, 
one obtains: 



dfSu ~ —1(5 



1 



dfSv = — -SuSv 



T 



l-lQTr 
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-^QTr 



Su 



Sv 



-QrAijriSvjSv ^ 



While the equation for longitudinal velocity increment is 
closed (if Q is known), it is apparent that the equation 
for transverse velocity increment cannot be closed in terms 
of Su, Sv and Q. As shown in [S], parts of the invariant 
Q may be expressed directly in terms of Su, namely Q = 
—SiP' I £'^+Q~ , where Q~ is composed of terms that may not 
be expressed in terms of Su or Sv. The argument above is 
useful, since it shows that within the approximation Q^ = 
(as was done in [5]), the added term from the linear ex- 
pansion of the pressure Hessian model is positive damping 
since (1 - 4Qrr/3) = {I + ASu^Tt/M'^) > 0, independent 
of Su. Hence, by setting Q^ — Q and neglecting the last 
term in the equation of the transverse velocity increment, 
the following system corresponds to the order RFD clo- 
sure applied to the ADV system (in three dimensions) : 




1 c 2 1 c 2 1 
— -Su^ + -Sv^ - — 

M £ T 
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(9) 



6. Conclusions and perspectives 

In this work we have explored several consequences of the 
RFD closure as applied to modeling pressure Hessian and 
viscous terms for the Lagrangian dynamics of the velocity 



n 



gradient tensor. Consistent with the observations of |6|7j , 
the model is shown to predict unphysical dynamics when 
attempting to reach high Reynolds numbers (decreasing 
the time-scale parameter t). In particular, for r < 10^^, 
results quickly deteriorate. Analysis of the time-correlation 
structure of the solution allowed us to conclude that an 
inconsistency develops, in which the assumed correlation 
time scale (r) becomes much smaller than the actual corre- 
lation time scale predicted by the dynamics of the system. 
Hence, at decreasing r, the assumptions underlying the clo- 
sure are not consistent with the simulated dynamics. How 
to "break" the observed longer-than-expected time corre- 
lations of the system is not clear and requires further study. 
Increasing the forcing strength was not a solution to the 
problem. In fact, a number of experiments were performed 
and in no case was it possible to obtain physically realis- 
tic predictions at high Reynolds number for r < 10"^. We 
must tentatively conclude that it may be impossible for 
a "single-shell" model in which only the dynamics at the 
smallest scales (largest velocity gradients) are computed 
dynamically to provide accurate predictions at arbitrarily 
high Reynolds numbers. As Re grows, it may be impossi- 
ble to describe the dynamics with only 8 independent de- 
grees of freedom, even though it is already remarkable that 
at moderate Re such a system is able to reproduce many 
features of turbulence [6|7j . This was the motivation for 
proposing a shell model [18] that included additional de- 
grees of freedom. 

With the aim at improving our understanding of the dy- 
namical effects of the recent fluid deformation closure, an 
expansion into powers of the small parameter (correlation 
time scale) was performed. It showed to be a good approx- 
imation, at least for r > 10~^. It allowed to show that, 
to first order, the RFD model for the pressure Hessian is 
a dissipative "damping" term for most of the tensor in- 
variant's time evolution, except for the second invariant of 
the strain-rate tensor in regions of high rotation. That is 
also the region in which direct comparisons with DNS in 
[S] showed significant inaccuracies of the closure. The effect 
of the linearized RFD closure for the pressure Hessian was 
shown to lead to a positive (dissipative) damping in the 
ADV equations © for longitudinal and transverse veloc- 
ity increments. For future efforts, it would be of interest to 
study the properties of this simple system of equations. In 
particular, one may want to develop stochastic forcing for 
the system in order to test whether it may yield realistic 
stationary statistics for longitudinal and transverse veloc- 
ity increments in turbulence. 



this article to Prof. K.R. Sreenivasan following the Sympo- 
sium on Fluid Science and Turbulence held on occasion of 
his eO*'' birthday. 

Appendix A. Tensorial random forcing term 

In this appendix, we describe in more detail the forcing 
F used in sections [3] and ID F is a random, traceless, ten- 
sorial noise that is restricted to comply with various tenso- 
rial symmetries to be consistent with those of the velocity 
gradient tensor. 

At first, we generate random deviates with uniform prob- 
ability distributions, and turn them into normal deviates 
by means of a simple transformation, specifying its Jaco- 
bian [IF . This is done for each of the 9 tensor components. 
Then, in order to obtain the white noise forcing used in 
some preliminary tests, and then as a comparison in figures 
[71- [TT|) . we adopt the procedure explained in appendix A of 
[5] . This procedure allows us to obtain the correct proper- 
ties for the trace of the tensor, and for the variances of the 
longitudinal and transverse components (the latter are the 
double of the former) . 

Next, we describe how to obtain the time-correlated 
noise. A time-uncorrelated tensorial noise {W) is used 
to obtain a finite-correlated [F) noise using an Ornstein- 
Uhlenbeck process [30] : 

dF^-6-^FAt + Te^'^dW , (A.l) 

where 9 and T are the desired correlation time and typical 
amplitude, respectively, for F . The correlation time will be 
set equal to r (tests with 9 = t/2, = t/5 and = t/10 were 
also presented, together with the white-noise case 9 = Q). 
As explained in section [31 the appropriate value of T to 
be prescribed is found empirically through figure [T] (and 
is corrected by introducing a multiplicative factor in sub- 
section [32] for various 9 values), and is such as to satisfy 
an energy-dissipation-based criterion. Note that the root- 
mean-square values of the transverse components are equal 
to JF, while a further factor l/\/2 is present for the longi- 
tudinal ones (e.g., F™^ ~ F™^ /^/i). Also, a single value 
of T must be used in (|A.ip for every component, because 
the correct ratio between the longitudinal and transverse 
amplitudes is guaranteed by the already-mentioned corre- 
sponding ratio for the white-noise components. 
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